Exploratory Spatial Analysis of Salmon Aquaculture

https://felixml.github.io/ECS530_h20/

Salmon Farm in Northern Norway

A large farm may contain as many fish as the entire population of wild salmon in Norway.

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tmap)
M <- tm_shape(Nor) + tm_polygons(, col= "grey")
P <- tm_shape(Akvaf) + tm_dots(col= "red" , size = 0.001 )
P2 <- tm_shape(Akvaf) + tm_dots(col= "red" , size = 0.01 )
M + P

With about 1,000 farms, the industry produces X1000 more than the natural carrying capacity of the ecosystem.

The biomass density of fish in the sea enables a parasite to thrive. The sea lice attaches to the skin of salmon and eats it, causing damages and infection. For the industry this is very costly. It also affects wild smolt migrating down the rivers into the sea.

The government has introduce several legislation which limits the growth of the industry.

library(raster)
st_bbox(Nor)
##      xmin      ymin      xmax      ymax 
##  3.904688 57.959026 31.161980 71.181252
rasta <- raster(Nor, st_bbox(Nor))
R1 <- rasterize(Akvaf,rasta, field = "LOK_NR", fun="count")
qtm(R1) + M

library(tmaptools)
Bomlo <- bb("Bomlo")

Bomlo2 <- tm_shape(Akvaf, bbox=Bomlo) +  tm_dots(col= "red" , size = 0.1 ) +
  M + tm_layout(bg.color="blue") + 
  tm_compass(type="8star", position = "left", color.dark = "white", color.light = "black", text.color="red" )

Bomlo2

library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:raster':
## 
##     extract
library(leaflet)
#Map Hardanger
tmap_leaflet(P2) %>% setView(lat=60.0, lng=5.8, zoom= 8)
#Map Bomlø
tmap_leaflet(P2) %>% setView(lat=59.8, lng=5.5, zoom= 9)

Nearest Neighbor

#Nearest Neighbor 1 : spdep
suppressPackageStartupMessages(library(spdep))
dnb <- dnearneigh(Akvaf, 0, 25, row.names=row.names(Akvaf))
dists <- nbdists(dnb, st_coordinates(Akvaf))
ozpo <- set.ZeroPolicyOption(TRUE)
lw <- nb2listw(dnb, glist=dists, style="B")
## Warning in nb2listw(dnb, glist = dists, style = "B"): zero sum general weights
#Nearest Neighbor 1 : nngeo
library(nngeo)
nn = st_nn(Akvaf, Akvaf, k = 3, progress = FALSE)
## lon-lat points
st_join(Akvaf, Akvaf, join = st_nn, k = 3, progress = FALSE)
## lon-lat points
## Simple feature collection with 3171 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 4.660783 ymin: 58.0243 xmax: 30.41533 ymax: 71.01778
## geographic CRS: WGS 84
## First 10 features:
##       LOK_NR.x LOK_NR.y                  geometry
## 157      12832    12832  POINT (8.26705 58.14398)
## 157.1    12832    14016  POINT (8.26705 58.14398)
## 157.2    12832    29957  POINT (8.26705 58.14398)
## 159      14016    14016   POINT (8.2553 58.14027)
## 159.1    14016    12832   POINT (8.2553 58.14027)
## 159.2    14016    29957   POINT (8.2553 58.14027)
## 163      29957    29957  POINT (7.06875 58.05583)
## 163.1    29957    30657  POINT (7.06875 58.05583)
## 163.2    29957    33577  POINT (7.06875 58.05583)
## 164      33577    33577 POINT (6.885683 58.06008)
l = st_connect(Akvaf, Akvaf, ids = nn, progress = FALSE)
plot(l, col = NA)  # For setting the extent
plot(st_geometry(Akvaf), col = "darkgrey", add = TRUE)
#plot(st_geometry(Akvaf), col = "red", add = TRUE)
plot(l, add = TRUE)

By Sea Distance

Buffers

#Buffers( What are the units here?)
buff25 <- st_buffer(Akvaf, 25)
## Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
## endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
## dist is assumed to be in decimal degrees (arc_degrees).
plot(buff25)
## Warning in st_is_longlat(x): bounding box has potentially an invalid value range
## for longlat data
## Warning in st_is_longlat(x): bounding box has potentially an invalid value range
## for longlat data

Lice :

https://www.barentswatch.no/nedlasting/fishhealth/lice

## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
summary(lus30$Lice)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.0000  0.0600  0.1208  0.1900  2.1200     488
nrow(lus30)
## [1] 1064
tapply(lus30$Lice, lus30$Fylke, summary)
## $Agder
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.00000 0.03000 0.06000 0.06333 0.09750 0.13000       7 
## 
## $`Møre og Romsdal`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.00000 0.00000 0.02000 0.06188 0.09000 0.36000      36 
## 
## $Nordland
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## 0.00000 0.00000 0.02000 0.09016 0.13000 0.78000     104 
## 
## $Rogaland
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.0225  0.1650  0.1947  0.3375  0.5000      35 
## 
## $`Trøndelag`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.0100  0.1000  0.1771  0.2500  2.1200      74 
## 
## $`Troms og Finnmark`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.0000  0.0200  0.0522  0.0500  0.4900     102 
## 
## $Vestland
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.0325  0.1300  0.1612  0.2700  0.4900     130

Vestland

#lm(lus30$Lice ~ lus30$Temp + I(lus30$Temp^2))
vest <- lus30[lus30$Fylke=="Vestland",]
nrow(vest)
## [1] 292
##    vars   n mean   sd median trimmed  mad min  max range skew kurtosis   se
## X1    1 162 0.16 0.14   0.13    0.15 0.15   0 0.49  0.49 0.62    -0.78 0.01

library(colorspace)
## 
## Attaching package: 'colorspace'
## The following object is masked from 'package:raster':
## 
##     RGB
pal <- heat_hcl(10)
mvest <- st_as_sf(vest, coords = c("Lon", "Lat"), crs = 4326)
plot(mvest["Lice"], nbreaks=7, breaks="fisher", col=pal) 
## Warning in plot.sf(mvest["Lice"], nbreaks = 7, breaks = "fisher", col = pal):
## col is not of length 1 or nrow(x): colors will be recycled; use pal to specify a
## color palette

How to Draw Density of Lice?

library(spatstat)
## Loading required package: spatstat.data
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## The following object is masked from 'package:raster':
## 
##     getData
## Loading required package: rpart
## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: spatstat version 1.64-1 is out of date by more than 6 months; we recommend upgrading to the latest version.
## 
## Attaching package: 'spatstat'
## The following object is masked from 'package:colorspace':
## 
##     coords
## The following objects are masked from 'package:psych':
## 
##     reflect, rescale
## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift
mvest <- st_as_sf(vest, coords = c("Lon", "Lat"), crs = 4326)
str(mvest)
## Classes 'sf' and 'data.frame':   292 obs. of  19 variables:
##  $ Uke                     : int  30 30 30 30 30 30 30 30 30 30 ...
##  $ Ã.r                     : int  2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
##  $ Lok_NB                  : int  30196 15196 10331 12067 12982 10298 11756 12086 13227 10319 ...
##  $ Lokalitetsnavn          : chr  "Ã…dnekvamme" "Aga Ã\230" "Ã…kre" "Aldalen" ...
##  $ Lice                    : num  NA 0.36 0.02 0.13 NA 0.18 NA 0.06 NA NA ...
##  $ Lus.i.bevegelige.stadier: num  NA 7.36 0.55 0.08 NA 0.32 NA 2.37 NA NA ...
##  $ Fastsittende.lus        : num  NA 1.73 0.56 0.18 NA 0.28 NA 0 NA NA ...
##  $ Brakklagt               : chr  "Ja" "Nei" "Nei" "Nei" ...
##  $ Telt                    : chr  "Nei" "Ja" "Ja" "Ja" ...
##  $ Kommunenummer           : int  4634 4613 4617 4624 4645 4615 4632 4617 4624 4602 ...
##  $ Kommune                 : chr  "MASFJORDEN" "BÃ\230MLO" "KVINNHERAD" "BJÃ\230RNAFJORDEN" ...
##  $ Fylkesnummer            : int  46 46 46 46 46 46 46 46 46 46 ...
##  $ Fylke                   : chr  "Vestland" "Vestland" "Vestland" "Vestland" ...
##  $ Lusegrense.uke          : chr  "0.5" "0.5" "0.5" "0.5" ...
##  $ Over.lusegrense.uke     : chr  "" "Nei" "Nei" "Nei" ...
##  $ Temp                    : num  NA 14.7 16.4 12.4 NA ...
##  $ ProduksjonsomrÃ.deId    : int  4 3 3 3 4 3 4 3 3 4 ...
##  $ ProduksjonsomrÃ.de      : chr  "Nordhordland til Stadt" "Karmøy til Sotra" "Karmøy til Sotra" "Karmøy til Sotra" ...
##  $ geometry                :sfc_POINT of length 292; first list element:  'XY' num  5.36 60.83
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:18] "Uke" "Ã.r" "Lok_NB" "Lokalitetsnavn" ...
#points$mark <- factor(mvest$Lice)
#points.ppp <- as.ppp(mvest)
#points.ppp$window <- as.owin(poly)

Interpolation : How to create variogram?

vest2 <- vest[vest$Telt=="JA",]
library(gstat)
## 
## Attaching package: 'gstat'
## The following object is masked from 'package:spatstat':
## 
##     idw
#hscat(Lice~Temp  , data=vest2)

g <- gstat(id="Lok_NB", formula=Lice ~ Temp, data=vest2)
str(g)
## List of 3
##  $ data :List of 1
##   ..$ Lok_NB:List of 15
##   .. ..$ formula      :Class 'formula'  language Lice ~ Temp
##   .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..$ data         :'data.frame': 0 obs. of  20 variables:
##   .. .. ..$ Uke                     : int(0) 
##   .. .. ..$ Ã.r                     : int(0) 
##   .. .. ..$ Lok_NB                  : int(0) 
##   .. .. ..$ Lokalitetsnavn          : chr(0) 
##   .. .. ..$ Lice                    : num(0) 
##   .. .. ..$ Lus.i.bevegelige.stadier: num(0) 
##   .. .. ..$ Fastsittende.lus        : num(0) 
##   .. .. ..$ Brakklagt               : chr(0) 
##   .. .. ..$ Telt                    : chr(0) 
##   .. .. ..$ Kommunenummer           : int(0) 
##   .. .. ..$ Kommune                 : chr(0) 
##   .. .. ..$ Fylkesnummer            : int(0) 
##   .. .. ..$ Fylke                   : chr(0) 
##   .. .. ..$ Lat                     : num(0) 
##   .. .. ..$ Lon                     : num(0) 
##   .. .. ..$ Lusegrense.uke          : chr(0) 
##   .. .. ..$ Over.lusegrense.uke     : chr(0) 
##   .. .. ..$ Temp                    : num(0) 
##   .. .. ..$ ProduksjonsomrÃ.deId    : int(0) 
##   .. .. ..$ ProduksjonsomrÃ.de      : chr(0) 
##   .. ..$ has.intercept: int 1
##   .. ..$ beta         : num(0) 
##   .. ..$ nmax         : num Inf
##   .. ..$ nmin         : num 0
##   .. ..$ omax         : num 0
##   .. ..$ maxdist      : num Inf
##   .. ..$ force        : logi FALSE
##   .. ..$ dummy        : logi FALSE
##   .. ..$ vfn          : int 1
##   .. ..$ weights      : NULL
##   .. ..$ degree       : num 0
##   .. ..$ vdist        : logi FALSE
##   .. ..$ lambda       : num 1
##  $ model: list()
##  $ call : language gstat(id = "Lok_NB", formula = Lice ~ Temp, data = vest2)
##  - attr(*, "class")= chr [1:2] "gstat" "list"
#evgm <- variogram(g)
#evgm <- variogram(g, cutoff=100, width=25)
#revgm <- variogram(g, cutoff=100, width=25, cressie=TRUE)
#cevgm <- variogram(g, cutoff=100, width=25, cloud=TRUE)

THANK YOU!!! :)